﻿
//说明：
//这一版 固定裁剪区域 固定模板区域 寻找坐标偏移量仅在一定范围内寻找（以模板在原图右下角坐标为基准画30*30方框 寻找对应图像中与模板最相似的部分的右下角点 进而得到坐标偏移量）
//可以达到每批均处理

#ifndef __CUDACC__
#define __CUDACC__
#endif // !__CUDACC__
#include "cuda_texture_types.h"

#include "cuda_runtime.h"
#include <thrust/extrema.h>
#include "device_launch_parameters.h"
#include<fstream>
#include <stdio.h>
#include<cuComplex.h>
#include<math.h>
#include <windows.h>
#include "opencv2/opencv.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc/types_c.h"
#include "opencv2/imgproc/imgproc_c.h"
#include "opencv2/highgui/highgui_c.h"
#include "kernel_stream_SAD_fixedarea.h"


#define CHECK(call){const cudaError_t error = call;if (error != cudaSuccess){printf("Error:%d\n", __LINE__);printf("code:%d,reason:%s\n", error, cudaGetErrorString(error));exit(1);}}
using namespace std;
using namespace cv;


//cudaEvent_t     start, stop;
//float elapsedTime;
//cudaEventCreate(&start);
//cudaEventCreate(&stop);
//cudaEventRecord(start, 0);
//cudaEventRecord(stop, 0);
//cudaEventSynchronize(stop);
//cudaEventElapsedTime(&elapsedTime, start, stop);
//cudaEventDestroy(start);
//cudaEventDestroy(stop);

texture<float, 1, cudaReadModeElementType> tex_line;
cudaTextureObject_t tex_gauss;


__global__ void GetData(float*data, float *complex, int*d_size, int*d_size_crop)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >=(d_size_crop[0]* d_size_crop[1]))
		return;
	int row_index = x % d_size_crop[0];//列内第几个
	int col_index = x / d_size_crop[0];//之前有多少列
	float Complex_temp = 0;
	int new_index_re = col_index * 2 * d_size[0] + row_index;
	int new_index_im = (col_index * 2 + 1)*d_size[0] + row_index;
	Complex_temp = log10f(sqrt(data[new_index_re]* data[new_index_re] + data[new_index_im]* data[new_index_im]));
	Complex_temp = (Complex_temp - 3) / 3;
	Complex_temp = Complex_temp < 0 ? 0 : Complex_temp;
	Complex_temp = Complex_temp > 1 ? 1 : Complex_temp;
	complex[x] = Complex_temp * 255;
}


__global__ void GetStru(float2*d_complex_crop, float* Bscan_ALL_GPU, int*d_size_crop)
{
	int ix = blockIdx.x*blockDim.x + threadIdx.x;
	if (ix >= (d_size_crop[0] * d_size_crop[1] * d_size_crop[2]))
		return;

	int frame_index = ix / (d_size_crop[0] * d_size_crop[1]);
	int row_index = ix % d_size_crop[0];
	int col_index = (ix - frame_index * (d_size_crop[0] * d_size_crop[1])) / d_size_crop[0];

	float *frameHead, *rowHead;
	frameHead = (float *)(Bscan_ALL_GPU + frame_index * d_size_crop[0] *  d_size_crop[1]);
	rowHead = (float*)(frameHead + col_index * d_size_crop[0]);
	rowHead[row_index] = log10f(cuCabsf(d_complex_crop[ix]) + 0.00000000000001);//加 存疑
}

__global__ void Gauss2D(float* d_gauss_data, int*d_size_crop, cudaTextureObject_t tex_gauss)//matlab中 3 3为中心 周围多余元素赋0
{
	int ix = blockIdx.x*blockDim.x + threadIdx.x;
	int iy = blockIdx.y*blockDim.y + threadIdx.y;
	int iz = blockIdx.z*blockDim.z + threadIdx.z;
	int index = (ix + iy * d_size_crop[0] + iz * (d_size_crop[0] * d_size_crop[1]));
	if (ix>= d_size_crop[0]|| iy >= d_size_crop[1]|| iz >= d_size_crop[2])
		return;

	/*float gauss_kernel[9] = { 0.0751,0.1238,0.0751,0.1238,0.2042,0.1238,0.0751,0.1238,0.0751 };
	d_gauss_data[index] = gauss_kernel[0] * tex3D<float>(tex_gauss, ix - 1, iy - 1, iz) + gauss_kernel[1] * tex3D<float>(tex_gauss, ix, iy - 1, iz) + gauss_kernel[2] * tex3D<float>(tex_gauss, ix + 1, iy - 1, iz) + gauss_kernel[3] * tex3D<float>(tex_gauss, ix - 1, iy, iz) + gauss_kernel[4] * tex3D<float>(tex_gauss, ix, iy, iz) +
		gauss_kernel[5] * tex3D<float>(tex_gauss, ix + 1, iy, iz) + gauss_kernel[6] * tex3D<float>(tex_gauss, ix - 1, iy + 1, iz) + gauss_kernel[7] * tex3D<float>(tex_gauss, ix, iy + 1, iz) + gauss_kernel[8] * tex3D<float>(tex_gauss, ix + 1, iy + 1, iz);*/

	float gauss_kernel[36] = { 0.0110073488022985,0.0181480501043652,0.0233025575973275,0.0233025575973275,0.0181480501043652,0.0110073488022985,0.0181480501043652,0.0299210762287985,0.0384194223724287,0.0384194223724287,0.0299210762287985,0.0181480501043652,
		0.0233025575973275,0.0384194223724287,0.0493315148206601,0.0493315148206601,0.0384194223724287,0.0233025575973275,0.0233025575973275,0.0384194223724287,0.0493315148206601,0.0493315148206601,0.0384194223724287,0.0233025575973275,
		0.0181480501043652,0.0299210762287985,0.0384194223724287,0.0384194223724287,0.0299210762287985,0.0181480501043652,0.0110073488022985,0.0181480501043652,0.0233025575973275,0.0233025575973275,0.0181480501043652,0.0110073488022985 };

	d_gauss_data[index] = gauss_kernel[0] * tex3D<float>(tex_gauss, ix - 2, iy - 2, iz) + gauss_kernel[1] * tex3D<float>(tex_gauss, ix - 1, iy - 2, iz) + gauss_kernel[2] * tex3D<float>(tex_gauss, ix, iy - 2, iz) + gauss_kernel[3] * tex3D<float>(tex_gauss, ix + 1, iy - 2, iz) + gauss_kernel[4] * tex3D<float>(tex_gauss, ix + 2, iy - 2, iz) + gauss_kernel[5] * tex3D<float>(tex_gauss, ix + 3, iy - 2, iz) +
		gauss_kernel[6] * tex3D<float>(tex_gauss, ix - 2, iy - 1, iz) + gauss_kernel[7] * tex3D<float>(tex_gauss, ix - 1, iy - 1, iz) + gauss_kernel[8] * tex3D<float>(tex_gauss, ix, iy - 1, iz) + gauss_kernel[9] * tex3D<float>(tex_gauss, ix + 1, iy - 1, iz) + gauss_kernel[10] * tex3D<float>(tex_gauss, ix + 2, iy - 1, iz) + gauss_kernel[11] * tex3D<float>(tex_gauss, ix + 3, iy - 1, iz) +
		gauss_kernel[12] * tex3D<float>(tex_gauss, ix - 2, iy, iz) + gauss_kernel[13] * tex3D<float>(tex_gauss, ix - 1, iy, iz) + gauss_kernel[14] * tex3D<float>(tex_gauss, ix, iy, iz) + gauss_kernel[15] * tex3D<float>(tex_gauss, ix + 1, iy, iz) + gauss_kernel[16] * tex3D<float>(tex_gauss, ix + 2, iy, iz) + gauss_kernel[17] * tex3D<float>(tex_gauss, ix + 3, iy, iz) +
		gauss_kernel[18] * tex3D<float>(tex_gauss, ix - 2, iy + 1, iz) + gauss_kernel[19] * tex3D<float>(tex_gauss, ix - 1, iy + 1, iz) + gauss_kernel[20] * tex3D<float>(tex_gauss, ix, iy + 1, iz) + gauss_kernel[21] * tex3D<float>(tex_gauss, ix + 1, iy + 1, iz) + gauss_kernel[22] * tex3D<float>(tex_gauss, ix + 2, iy + 1, iz) + gauss_kernel[23] * tex3D<float>(tex_gauss, ix + 3, iy + 1, iz) +
		gauss_kernel[24] * tex3D<float>(tex_gauss, ix - 2, iy + 2, iz) + gauss_kernel[25] * tex3D<float>(tex_gauss, ix - 1, iy + 2, iz) + gauss_kernel[26] * tex3D<float>(tex_gauss, ix, iy + 2, iz) + gauss_kernel[27] * tex3D<float>(tex_gauss, ix + 1, iy + 2, iz) + gauss_kernel[28] * tex3D<float>(tex_gauss, ix + 2, iy + 2, iz) + gauss_kernel[29] * tex3D<float>(tex_gauss, ix + 3, iy + 2, iz) +
		gauss_kernel[30] * tex3D<float>(tex_gauss, ix - 2, iy + 3, iz) + gauss_kernel[31] * tex3D<float>(tex_gauss, ix - 1, iy + 3, iz) + gauss_kernel[32] * tex3D<float>(tex_gauss, ix, iy + 3, iz) + gauss_kernel[33] * tex3D<float>(tex_gauss, ix + 1, iy + 3, iz) + gauss_kernel[34] * tex3D<float>(tex_gauss, ix + 2, iy + 3, iz) + gauss_kernel[35] * tex3D<float>(tex_gauss, ix + 3, iy + 3, iz);
}

__global__ void CropRefData(float* d_gauss_data, float* d_ref_data, int*d_location, int*d_size_Ref, int*d_size_crop)
{
	int ix = blockIdx.x*blockDim.x + threadIdx.x;
	int z_num = d_size_crop[0];
	int x_num = d_size_crop[1];
	int y_num = d_size_crop[2];

	if (ix >= d_size_Ref[0] * d_size_Ref[1] * d_size_Ref[2])
		return;

	int frame_index = ix / (d_size_Ref[0] * d_size_Ref[1]);
	int row_index = ix % d_size_Ref[0];
	int col_index = (ix % (d_size_Ref[0] * d_size_Ref[1])) / d_size_Ref[0];

	d_ref_data[ix] = d_gauss_data[frame_index * 3 * z_num*x_num + (col_index + d_location[0])*z_num + row_index + d_location[2]];
}

__global__ void GetSAD(float*d_Frame_ROI_Ref, float*d_FrameRef, float*d_SAD, int *d_size_Frame_ROI_Ref, int *d_size_FrameRef,int*d_size_traverse, int *d_location)
{
	int ix = blockDim.x*blockIdx.x + threadIdx.x;
	int iy = blockDim.y*blockIdx.y + threadIdx.y;

	float sum = 0;
	int size_Frame_ROI_Ref_cols = d_size_Frame_ROI_Ref[1];
	int size_Frame_ROI_Ref_rows = d_size_Frame_ROI_Ref[0];
	int size_FrameRef_cols = d_size_FrameRef[1];
	int size_FrameRef_rows = d_size_FrameRef[0];

	int size_SAD_rows = d_size_traverse[0];
	if (ix >= d_size_traverse[0] || iy >= d_size_traverse[1])
		return;

	int index_x = d_size_Frame_ROI_Ref[0] - 1 + d_location[2] - d_size_traverse[0] / 2 + ix;//模板右下角坐标在遍历区域的中心
	int index_y = d_size_Frame_ROI_Ref[1] - 1 + d_location[0] - d_size_traverse[1] / 2 + iy;

	for (int i = 0; i < size_Frame_ROI_Ref_cols; i++)
	{
		for (int j = 0; j < size_Frame_ROI_Ref_rows; j++)
		{
			sum += fabsf(d_Frame_ROI_Ref[(size_Frame_ROI_Ref_cols - i - 1)*size_Frame_ROI_Ref_rows + (size_Frame_ROI_Ref_rows - j - 1)] - d_FrameRef[(index_y - i)*size_FrameRef_rows + index_x - j ]);
		}
	}

	d_SAD[iy*size_SAD_rows + ix] = sum;
}

__global__ void GetMin(float*d_SAD, int2* d_coord, int* d_size_traverse, int *d_size_Frame_ROI_Ref, int *d_location)
{
	int coord = 0;
	float min = d_SAD[0];
	int traverse_size_all = d_size_traverse[0] * d_size_traverse[1];
	for (int i = 1; i < traverse_size_all; i++)
	{
		float SAD = d_SAD[i];
		if (SAD < min)
		{
			min = SAD;
			coord = i;
		}
	}
	d_coord[0].x = (int)(coord % d_size_traverse[0] - d_size_traverse[0] / 2);//- (d_size_Frame_ROI_Ref[0]-1+ d_location[2])
	d_coord[0].y = (int)(coord / d_size_traverse[0] - d_size_traverse[1] / 2);//-( d_size_Frame_ROI_Ref[1]-1+9)
}

__global__ void GetAligned(float2*d_complex_crop, float2*d_complex_crop_aligned, int2* d_coord,int*d_size_crop)
{
	int ix = blockDim.x*blockIdx.x + threadIdx.x;
	int iy = blockDim.y*blockIdx.y + threadIdx.y;
	if (ix >= d_size_crop[0] || iy >= d_size_crop[1])
		return;

	int x_border_1 = 0;
	int x_border_2 = d_size_crop[0] - 1;
	int y_border_1 = 0;
	int y_border_2 = d_size_crop[1] - 1;
	int x_offset = d_coord[0].x;
	int y_offset = d_coord[0].y;
	
	if (x_offset < 0)
		x_border_1 = -x_offset;//小于x_border_1的均赋0
	else
		x_border_2 = d_size_crop[0] - 1 - x_offset;//大于x_border_2的均赋0

	if (y_offset < 0)
		y_border_1 = -y_offset;//小于x_border_1的均赋0
	else
		y_border_2 = d_size_crop[1] - 1 - y_offset;//大于x_border_2的均赋0
	
	int index = ix + iy * d_size_crop[0];
	int index_origin = ix + x_offset + (iy+y_offset) * d_size_crop[0];

	if (ix< x_border_1 || ix>x_border_2 || iy< y_border_1 || iy>y_border_2)
	{
		d_complex_crop_aligned[index].x = 0;
		d_complex_crop_aligned[index].y = 0;
	}
	else
	{
		d_complex_crop_aligned[index].x = d_complex_crop[index_origin].x;
		d_complex_crop_aligned[index].y = d_complex_crop[index_origin].y;
	}

}

__global__ void preImfilter_C1(cuFloatComplex*complex, cuFloatComplex*preImfilter, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3 * 2)
		return;
	int frame_index = x / d_size[0];
	int new_frame_index = frame_index / 2 * 3 + frame_index % 2 + 1;
	int index = new_frame_index * d_size[0] + x % d_size[0];//两帧映射回原来三帧的坐标

	int size_last_frame = d_size[0] * (d_size[1] - 1);
	int size_first_frame = d_size[0];
	int shift_index = 0;//complex_shift
	if (index < size_first_frame)
		shift_index = index + size_last_frame;
	else
		shift_index = index - d_size[0];
	preImfilter[x] = cuCmulf(complex[index], cuConjf(complex[shift_index]));//注意原始坐标对应 complex[shift_index]即complex中哪项对应complex_shift[x]

}

__global__ void preImfilter_C2(cuFloatComplex*complex, float*preImfilter, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3 * 2)
		return;
	int frame_index = x / d_size[0];
	int new_frame_index = frame_index / 2 * 3 + frame_index % 2 + 1;
	int index = new_frame_index * d_size[0] + x % d_size[0];//两帧映射回原来三帧的坐标
	float abs = cuCabsf(complex[index]);
	preImfilter[x] = abs * abs;
}

__global__ void preImfilter_C3(cuFloatComplex*complex, float*preImfilter, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3 * 2)
		return;

	int frame_index = x / d_size[0];
	int new_frame_index = frame_index / 2 * 3 + frame_index % 2 + 1;
	int index = new_frame_index * d_size[0] + x % d_size[0];//两帧映射回原来三帧的坐标

	int size_last_frame = d_size[0] * (d_size[1] - 1);
	int size_first_frame = d_size[0];
	int shift_index = 0;//complex_shift
	if (index < size_first_frame)
		shift_index = index + size_last_frame;
	else
		shift_index = index - d_size[0];
	float abs = cuCabsf(complex[shift_index]);
	preImfilter[x] = abs * abs;
}

__global__ void Imfilter_hC(float*preImfilter, float*c2, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3 * 2)
		return;
	int col_index = x % d_size[2];
	float a, b, c, d, e;

	if (col_index > 1 && col_index < d_size[2] - 3)
	{
		a = preImfilter[x - 2];
		b = preImfilter[x - 1];
		c = preImfilter[x + 1];
		d = preImfilter[x + 2];
		e = preImfilter[x + 3];
	}
	else if (col_index == 1)
	{
		a = preImfilter[x - 1];
		b = preImfilter[x - 1];
		c = preImfilter[x + 1];
		d = preImfilter[x + 2];
		e = preImfilter[x + 3];
	}
	else if (col_index == d_size[2] - 3)
	{
		a = preImfilter[x - 2];
		b = preImfilter[x - 1];
		c = preImfilter[x + 1];
		d = preImfilter[x + 2];
		e = preImfilter[x + 2];
	}
	else if (col_index == d_size[2] - 2)
	{
		a = preImfilter[x - 2];
		b = preImfilter[x - 1];
		c = preImfilter[x + 1];
		d = preImfilter[x + 1];
		e = preImfilter[x + 1];
	}
	else if (col_index == d_size[2] - 1)
	{
		a = preImfilter[x - 2];
		b = preImfilter[x - 1];
		c = preImfilter[x];
		d = preImfilter[x];
		e = preImfilter[x];
	}
	else if (col_index == 0)
	{
		a = preImfilter[x];
		b = preImfilter[x];
		c = preImfilter[x + 1];
		d = preImfilter[x + 2];
		e = preImfilter[x + 3];
	}

	//由于switch语句是查找表 需要在编译阶段确定case 不能传入 暂时不用switch

	c2[x] = sqrt(a + b + preImfilter[x] + c + d + e);
	return;
}

__global__ void Imfilter_h(cuFloatComplex*complex, cuFloatComplex*c1, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3 * 2)
		return;
	int col_index = x % d_size[2];
	cuFloatComplex a, b, c, d, e;

	if (col_index > 1 && col_index < d_size[2] - 3)
	{
		a = complex[x - 2];
		b = complex[x - 1];
		c = complex[x + 1];
		d = complex[x + 2];
		e = complex[x + 3];
	}
	else if (col_index == 1)
	{
		a = complex[x - 1];
		b = complex[x - 1];
		c = complex[x + 1];
		d = complex[x + 2];
		e = complex[x + 3];
	}
	else if (col_index == d_size[2] - 3)
	{
		a = complex[x - 2];
		b = complex[x - 1];
		c = complex[x + 1];
		d = complex[x + 2];
		e = complex[x + 2];
	}
	else if (col_index == d_size[2] - 2)
	{
		a = complex[x - 2];
		b = complex[x - 1];
		c = complex[x + 1];
		d = complex[x + 1];
		e = complex[x + 1];
	}
	else if (col_index == d_size[2] - 1)
	{
		a = complex[x - 2];
		b = complex[x - 1];
		c = complex[x];
		d = complex[x];
		e = complex[x];
	}
	else if (col_index == 0)
	{
		a = complex[x];
		b = complex[x];
		c = complex[x + 1];
		d = complex[x + 2];
		e = complex[x + 3];
	}

	c1[x].x = a.x + b.x + complex[x].x + c.x + d.x + e.x;
	c1[x].y = a.y + b.y + complex[x].y + c.y + d.y + e.y;
	return;
}

__global__ void GetC(cuFloatComplex*c1, float*c2, float*c3, float*c, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3 * 2)
		return;

	float a = c2[x] * c3[x];
	cuFloatComplex b = c1[x];
	b.x = b.x / a;
	b.y = b.y / a;
	c[x] = cuCabsf(b);
}

__global__ void GetD(float*c, float*d, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3)
		return;
	int frame_index = x / d_size[0];
	int index_a = frame_index * 2 * d_size[0] + x % d_size[0];
	int index_b = (frame_index * 2 + 1) * d_size[0] + x % d_size[0];
	d[x] = 1 - 0.5*(c[index_a] + c[index_b]);
}

__global__ void GetStruAndiSNR(cuFloatComplex*complex, float*iSNR, float* Stru, float s2, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3)
		return;
	int frame_index = x / d_size[0];
	int index_a = frame_index * 3 * d_size[0] + x % d_size[0];
	int index_b = (frame_index * 3 + 1) * d_size[0] + x % d_size[0];
	int index_c = (frame_index * 3 + 2) * d_size[0] + x % d_size[0];
	float a = cuCabsf(complex[index_a]);
	float b = cuCabsf(complex[index_b]);
	float c = cuCabsf(complex[index_c]);

	Stru[x] = log10f((a + b + c) / 3);
	iSNR[x] = s2 * 3 / (a*a + b * b + c * c);
}

__global__ void GetidDecor(float*d, float*iSNR, float* idDecor, float isnr_Coefficient, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3)
		return;
	if (d[x] > iSNR[x] * isnr_Coefficient)
		idDecor[x] = d[x];
	else
		idDecor[x] = 0;
}

__global__ void GetProjection(float* idDecor, float* projection, int*d_size, int *d_range)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= (d_size[1] * (d_size[0] / d_size[2]) / 3))
		return;
	int frame_num = x / (d_size[0] / d_size[2]);
	int col_index = x % (d_size[0] / d_size[2]);
	int index = frame_num * d_size[0] + col_index * d_size[2];

	int range_high = d_range[0];
	if (range_high >= d_size[2])
		range_high = 0;
	int range_low = d_range[1];
	if (range_low >= d_size[2])
		range_low = d_size[2]-1;

	float max = idDecor[index + d_range[0]];
	for (int i = range_high; i < range_low; i++)
	{
		float a = idDecor[index + i];
		if (a> max)
			max = a;
	}
	projection[x] = max;
}

__global__ void GaussRow(float* idDecor, float* result, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3)
		return;
	int frame_index = x % d_size[0];//帧内坐标
	int frame_num = x / d_size[0];//之前有几帧
	int row_index = frame_index % d_size[2];//第几行
	int col_index = frame_index / d_size[2];//之前有几列
	int left_index = frame_num * d_size[0] + (col_index - 1)*d_size[2] + row_index;
	int right_index = frame_num * d_size[0] + (col_index + 1)*d_size[2] + row_index;
	if (col_index > 0 && col_index < ((d_size[0] / d_size[2]) - 1))
		result[x] = idDecor[x] * 0.7870 + idDecor[left_index] * 0.1065 + idDecor[right_index] * 0.1065;
	else if (col_index == 0)
		result[x] = idDecor[x] * 0.7870 + idDecor[x] * 0.1065 + idDecor[right_index] * 0.1065;
	else if (col_index == (d_size[0] / d_size[2] - 1))
		result[x] = idDecor[x] * 0.7870 + idDecor[left_index] * 0.1065 + idDecor[x] * 0.1065;
}

__global__ void GaussCol(float* gaussRow, float* result, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3)
		return;
	int col_index = x % d_size[2];
	if (col_index > 3 && col_index < (d_size[2] - 4))

		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x - 1] + gaussRow[x + 1]) + 0.0219*(gaussRow[x - 2] + gaussRow[x + 2]) + (4.4074e-4)*(gaussRow[x - 3] + gaussRow[x + 3]) + (1.8584e-6)*(gaussRow[x - 4] + gaussRow[x + 4]);
	else if (col_index == 0)
		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x] + gaussRow[x + 1]) + 0.0219*(gaussRow[x] + gaussRow[x + 2]) + (4.4074e-4)*(gaussRow[x] + gaussRow[x + 3]) + (1.8584e-6)*(gaussRow[x] + gaussRow[x + 4]);
	else if (col_index == 1)
		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x - 1] + gaussRow[x + 1]) + 0.0219*(gaussRow[x - 1] + gaussRow[x + 2]) + (4.4074e-4)*(gaussRow[x - 1] + gaussRow[x + 3]) + (1.8584e-6)*(gaussRow[x - 1] + gaussRow[x + 4]);
	else if (col_index == 2)
		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x - 1] + gaussRow[x + 1]) + 0.0219*(gaussRow[x - 2] + gaussRow[x + 2]) + (4.4074e-4)*(gaussRow[x - 2] + gaussRow[x + 3]) + (1.8584e-6)*(gaussRow[x - 2] + gaussRow[x + 4]);
	else if (col_index == 3)
		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x - 1] + gaussRow[x + 1]) + 0.0219*(gaussRow[x - 2] + gaussRow[x + 2]) + (4.4074e-4)*(gaussRow[x - 3] + gaussRow[x + 3]) + (1.8584e-6)*(gaussRow[x - 3] + gaussRow[x + 4]);
	else if (col_index == (d_size[2] - 4))
		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x - 1] + gaussRow[x + 1]) + 0.0219*(gaussRow[x - 2] + gaussRow[x + 2]) + (4.4074e-4)*(gaussRow[x - 3] + gaussRow[x + 3]) + (1.8584e-6)*(gaussRow[x - 4] + gaussRow[x + 3]);
	else if (col_index == (d_size[2] - 3))
		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x - 1] + gaussRow[x + 1]) + 0.0219*(gaussRow[x - 2] + gaussRow[x + 2]) + (4.4074e-4)*(gaussRow[x - 3] + gaussRow[x + 2]) + (1.8584e-6)*(gaussRow[x - 4] + gaussRow[x + 2]);
	else if (col_index == (d_size[2] - 2))
		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x - 1] + gaussRow[x + 1]) + 0.0219*(gaussRow[x - 2] + gaussRow[x + 1]) + (4.4074e-4)*(gaussRow[x - 3] + gaussRow[x + 1]) + (1.8584e-6)*(gaussRow[x - 4] + gaussRow[x + 1]);
	else if (col_index == (d_size[2] - 1))
		result[x] = 0.4987*gaussRow[x] + 0.2283*(gaussRow[x - 1] + gaussRow[x]) + 0.0219*(gaussRow[x - 2] + gaussRow[x]) + (4.4074e-4)*(gaussRow[x - 3] + gaussRow[x]) + (1.8584e-6)*(gaussRow[x - 4] + gaussRow[x]);
}

__global__ void GaussSlc(float* gaussCol, float* result, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3)
		return;

	int frame_num = x / d_size[0];
	if (frame_num > 0 && frame_num < (d_size[1] / 3 - 1))//注意血流图帧数缩减为原来三分之一
		result[x] = 0.7870*gaussCol[x] + 0.1065*(gaussCol[x - d_size[0]] + gaussCol[x + d_size[0]]);

	else if (frame_num == 0)
		result[x] = 0.7870*gaussCol[x] + 0.1065*(gaussCol[x] + gaussCol[x + d_size[0]]);
	else if (frame_num == (d_size[1] / 3 - 1))
		result[x] = 0.7870*gaussCol[x] + 0.1065*(gaussCol[x - d_size[0]] + gaussCol[x]);
}

__global__ void medfilt3(float* gaussSlc, float* result, int*d_size)
{
	int x = blockIdx.x*blockDim.x + threadIdx.x;
	if (x >= d_size[0] * d_size[1] / 3)
		return;
	int col_index = x % d_size[2];
	float med;
	if (col_index > 0 && col_index < (d_size[2] - 1))
	{
		if ((gaussSlc[x] <= gaussSlc[x + 1] && gaussSlc[x] >= gaussSlc[x - 1]) || (gaussSlc[x] <= gaussSlc[x - 1] && gaussSlc[x] >= gaussSlc[x + 1]))
			med = gaussSlc[x];
		else if ((gaussSlc[x - 1] <= gaussSlc[x] && gaussSlc[x - 1] >= gaussSlc[x + 1]) || (gaussSlc[x - 1] <= gaussSlc[x + 1] && gaussSlc[x - 1] >= gaussSlc[x]))
			med = gaussSlc[x - 1];
		else
			med = gaussSlc[x + 1];
	}
	else
		med = gaussSlc[x];//采用'replicate' 与源程序'symmetric'不同
	result[x] = med;
}


__global__ void Avefilter(float* stru, float* stru_ave, int* d_size)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	int z = blockIdx.z * blockDim.z + threadIdx.z;

	int nx = d_size[0];
	int ny = d_size[1];
	int nz = d_size[2];

	if (x > nx - 1 || y > ny - 1 || z > nz - 1)
		return;

	int index_ave = z * nx * ny + y * nx + x;

	if (x<2 || x >nx - 3)
		stru_ave[index_ave] = stru[index_ave];
	else if (y<2 || y > ny - 3)
	{
		float sum = 0;
		for (int j = -2; j < 3; j++)
		{
			int index = z * nx * ny + y * nx + (x + j);
			sum += stru[index];
		}
		stru_ave[index_ave] = sum / 5.0f;
	}
	else
	{
		float sum = 0;
		for (int i = -2; i < 3; i++)
		{
			for (int j = -2; j < 3; j++)
			{
				int index = z * nx * ny + (y + i) * nx + (x + j);
				sum += stru[index];
			}
		}
		stru_ave[index_ave] = sum / 25.0f;
	}
}


__global__ void Getgradient(float* stru_ave, float* gradient, int* d_size)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;
	int y = blockIdx.y * blockDim.y + threadIdx.y;
	int z = blockIdx.z * blockDim.z + threadIdx.z;

	int nx = d_size[0];
	int ny = d_size[1];
	int nz = d_size[2];

	if (x > nx - 1 || y > ny - 1 || z > nz - 1)
		return;

	int index_a = z * nx * ny + y * nx + x;

	if (x<1 || x >nx - 2)
		gradient[index_a] = 0;
	else
	{
		gradient[index_a] = 0.5 * (stru_ave[index_a + 1] - stru_ave[index_a - 1]);
	}
}


__global__ void GetLayerTemp(int* d_gradient_location, float* gradient, int* d_size, int start_location, int end_location)
{
	int x = blockIdx.x * blockDim.x + threadIdx.x;

	int nx = d_size[0];
	int ny = d_size[1];
	int nz = d_size[2];

	if (x > ny * nz - 1)
		return;

	int index_z = x / ny;
	int start_index = x * nx + start_location;
	int end_index = x * nx + end_location;

	float max_a = 0;
	int location_a = start_location;

	for (int i = start_index; i < end_index; i++)
	{
		if (gradient[i] > max_a)
		{
			max_a = gradient[i];
			location_a = i - x * nx;
		}
	}

	d_gradient_location[x] = location_a;

}



void DataProcess(void* data[12], int* h_size, float* h_stru)
{
	int z_num = h_size[0];// A_scan/2
	int x_num = h_size[1];// B_scan
	int y_num = h_size[2];// 12
	
	const int size = h_size[0] * h_size[1] * h_size[2];
	int *d_size = nullptr;
	CHECK(cudaMalloc((void**)&d_size, 3 * sizeof(int)));
	CHECK(cudaMemcpy(d_size, h_size, 3 * sizeof(int), cudaMemcpyHostToDevice));


	//从alazar中读入数据 数据裁剪 暂时取消裁剪 使用固定大小300
	int size_crop[3] = { 512 ,x_num ,y_num };
	int crop_frame_size = size_crop[0]* size_crop[1];
	int *d_size_crop = nullptr;
	CHECK(cudaMalloc((void**)&d_size_crop, 3 * sizeof(int)));
	CHECK(cudaMemcpy(d_size_crop, size_crop, 3 * sizeof(int), cudaMemcpyHostToDevice));

	int size_crop_all = size_crop[0] * size_crop[1] * size_crop[2];
	float *d_complex_crop = nullptr;
	CHECK(cudaMalloc((void**)&d_complex_crop, size_crop_all * sizeof(float)));

	// ADD

	const dim3 blockSizeGetData(64);
	const dim3 gridSizeGetData((crop_frame_size + blockSizeGetData.x - 1) / blockSizeGetData.x);
	
	for (int i = 0; i < y_num; i++)
	{
		float* b = (float*)data[i];
		GetData << < gridSizeGetData, blockSizeGetData >> > (b, d_complex_crop+ i*crop_frame_size, d_size, d_size_crop);
	}
	

	CHECK(cudaMemcpy(h_stru, d_complex_crop, size_crop_all * sizeof(float), cudaMemcpyDeviceToHost));

	//结构图以及高斯滤波  4ms
	//const dim3 blockSizeCro(32);
	//const dim3 gridSizeCro((size_crop_all + blockSizeCro.x - 1) / blockSizeCro.x);

	//float* Bscan_ALL_GPU = nullptr;
	//CHECK(cudaMalloc((void**)&Bscan_ALL_GPU, size_crop_all * sizeof(int)));
	//GetStru << < gridSizeCro, blockSizeCro >> > (d_complex_crop, Bscan_ALL_GPU, d_size_crop);

	

	CHECK(cudaFree(d_size));
	CHECK(cudaFree(d_size_crop));
	CHECK(cudaFree(d_complex_crop));
	//CHECK(cudaFree(Bscan_ALL_GPU));

}



